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Abstract 


We study a diffusive energy -balance climate model, 
governed by a nonlinear parabolic partial differential 
equation. Three positive steady-state solutions of this 
equation are found? they correspond to three possible 
climates of our planet: an interglacial (nearly identical 

to the present climate) , a glacial, and a completely ice- 
covered earth. We consider also models similar to the main 
one studied, and determine the number of their steady states. 
All the models have albedo continuously varying with latitude 
and temperature, and entirely diffusive horizontal heat 
transfer. The diffusion is taken to be nonlinear as well as 
linear. 

We investigate the stability under small perturbations 
of the main model’s climates. A stability criterion is 
derived, and its application shows that the "present climate" 
and the "deep freeze" are stable, whereas the model's glacial 
is unstable. A variational principle is introduced to confirm 
the results of this stability analysis. 

We examine the dependence of the number of steady states 
and of their stability on the average solar radiation. The 
main result is that for a sufficient decrease in solar radia- 
tion (about 2 percent) the glacial and interglacial solutions 
disappear, leaving the ice-covered earth as the only possible 
climate . 
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1. Introduction 


The concept of a climate is one of those abstractions 
which appears to be self-evident to the layman, but is by 
no means well defined scientifically. The intuitive idea 
of a climate has two aspects: 

(a) the most important features of atmospheric phenomena; 

(b) the average behavior of these phenomena over a suitably 
long time interval and over sufficiently large areas. 

The difficulties start when one tries to give a precise 

meaning to the key words "most important", "suitably long" 
and "suitably large". We start with "suitably long"; 
clearly, a year is an absolute lower bound for a reasonable 
averaging time interval, since daily and seasonal variations 
should be excluded. To decide over how much longer than a 
year the averaging should be performed, one has to look at 
the record. There are three kinds of records: instrumental, 

the length of which is of the order of hundreds of years, 
historical, of the order of thousands of years, and geological, 
of the order of hundreds of thousands of years. These 
records show that features of the atmosphere change on all 
the time scales represented in them (e. g. , Robinson, 1971), 

Thus it would appear at first that it is not possible to 
distinguish between "fast" variations in yearly averages, 
which should be averaged out when defining a climate, and 
"slow" variations, which should be considered as "changes 
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of climate". Still, the geological record seems to indicate 
that the transitions between considerably colder periods (ice 
ages or glacials) and warmer periods ("normal" climates or 
interglacials) occurred over time spans about ten times 
shorter than the duration of the relatively steady cold or 
warm weather respectively. This suggests what we . shall adopt 
here as our operative definition of climate., viz., the preva- 
lence of either warm weather (as we experience it today) or 
of cold weather (to mean a difference of the order of ten 
degrees centigrade in yearly average below the one recorded 
in the present) . 

We. turn now to the question of which features of 
atmospheric phenomena are "most important" . Certainly tempera- 
ture is one of them, not only because its changes left deep 
traces in the geological record (glaciations in temperate 
zones, pluviations in the tropics — SM1C, 1971), but also 
because it affects all conditions of life and because it is 
directly linked to the major thermodynamic and dynamic 
processes in the atmosphere which determine climate and 
its changes. Also, humidity, wind direction and intensity, 
cloud amount, precipitation, all play a major 
role in determining what is perceived as weather and hence 
should be time-averaged (and, possibly, space-averaged) 
into climate. Moreover, it is not only the averages of 
these quantities, but their day-to-night and season-to-season 
contrast that enters our intuitive concept of a climate. 
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Thus, at least their variance should be included in a more 
complete mathematical model for climatology. 

In this article we shall treat a very simple model , 
based on the work of Sellers (1969) and of Schneider 
and Gal-Chen (1973); we hope that the results will in 
themselves be of some significance for climate theory, as 
well as providing Insights for devising and analyzing more 
complex models. 

In Section 2 the model to be studied is described; 
the physical principles on which it is based, as well as 
the empirical data it uses are discussed. 

In Section 3 we discuss the work of different authors 
on similar models; the similarities and differences between 
their results are pointed out and the issues arising from 
these results are outlined. . 

In Section 4 we compute numerically the model's 
steady-state solutions of physical interest, i.e., those 
yielding positive absolute temperatures. Three such 
solutions, corresponding to three distinct climates of our 
planet, are obtained; one corresponds to the current climate, 
the second to an ice age * the third to a completely ice- 
covered earth. In this section we also explore the effect 
of certain modifications in the model on the number of 
steady-state solutions. 

In Section 5 the notion of stability for the solutions 
obtained in Section 4 is defined precisely; it is investigated 
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using a combination of analytical and numerical techniques. 
The results are that the present climate and the ice-covered 
earth are stable, whereas the ice age of the model is 
unstable . 

In Section 6 the effect of changes in the solar radia- 
tion on t he n umb e r and stab ilit y of s te ady-s tate_s.pl u t ion s _ 
is studied. The main result is that for a sufficient 
decrease in the solar radiation (about 2%) , the glacial and 
the interglacial solutions disappear, leaving the ice- 
covered earth as the only possible climate. 
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2. The Model 


The model chosen for study is a zonally-and-vertically 
averaged energy -balance climate model. This means that 
quantities in the model are averaged over longitude and 
height, leaving colatitude 4> as the only space variable. 
The term energy balance means that the model is essentially 
based on the energy equation of fluid dynamics and has sea- 
level temperature u as the only dependent variable. The 
equation governing the model is 

(la) C ( 4? ) u t = R i (<f>,u) - R q ( cf> , u) + D (<j>,u,u^,u^) ; 

C is the heat capacity of the atmosphere, land and water 
masses; is the heat absorbed from incoming radiation, 

(lb) - G(4>) [1 - a (4>,u) ] , 

v. 

where Q is high-frequency solar radiation and a is the 
reflectivity (albedo) of the land and sea surface; Rq is 
the heat lost in outgoing low-frequency planetary re- 
radiation reaching outer space, 

4 

(lc) R q = c ( <f) , u) cru ; 

and D describes the redistribution of heat on the surface 
of the planet by conduction and convection, 

< ld) D ■ sdr? V sin k( +' u) i u * • 

The coefficients and forcing terms in this model 
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represent yearly averages of the corresponding quantities and 
therefore do not depend explicitly on the time t. Averaging 
the daily and seasonal variations seems justified, since the 
time scales in which we are interested in our investigation 
are of the order of hundreds and thousands of years . The lack 
of explicit dependence on time has the advantage that the 
model admits, as we shall see, steady-state solutions, u^ = 0, 
which we define to be its climates . The purpose of this work 
is to study analytically and numerically the number of these 
climates and their stability under perturbations. 

The first model of this type, in finite-difference form 
and without time dependence, was developed by Sellers (1969) . 
The differential formulation is due to Faegre (1972); time 
dependence was introduced by Dwyer and Petersen (1973) and, 
independently, by Schneider and Gal-Chen (1973) . Dwyer and 
Petersen also gave the outline of a systematic derivation of 
(1) from the energy equation of fluid dynamics, mentioning the 
main assumptions involved. An even simpler model has been 
proposed by Budyko (1969): in it the diffusion term D is 
replaced by a nondif ferentiated, linear term in u, and the 
albedo is a simple step function of u only; this model was 
also discussed very thoroughly by Leith (1974), and by Held 
and Suarez (1974). 

One of the main features of the model (la-d) is the form 
of the albedo, 

(2a) a = {b ( <f> ) - CjJu^ + (u-c 2 z (<j>) -u ) _] > c , 

where the meaning of the subscripts ( ) _ and { } is 

“ c 
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given for a generic quantity h by 
(2b) h_ = min {h,0} 


and 


(2c) 


h 


c 


f 0.25, h £ 0.25, 
• h , 0.25 < h < 0. 85, 
*■0.85,, 0.85 < h 


The subscript c stands for cutoff? the cutoff given by (2c) 
embodies the observed minimum and maximum values of surface 
albedo. 

Snow and ice have higher reflectivity than bare ground 
or water; since in regions of lower yearly average tempera- 
ture the snow and ice cover persists for a longer fraction 
of the year, at lower temperature the yearly average albedo 
is higher; this is expressed in the monotonically decreasing 
dependence of a on u. Further, the plausible assumption is 
made that, above a certain yearly average temperature u m , 
no snow or ice will be present at any time of the year; 
therefore a is independent of u for u-c 2 z > , as seen 

from (2a) and (2b). The term c 2 z(<J)) gives the difference 
between sea-level temperature u and ground temperature, u-c 2 z. 

A serious drawback of the model is that it does not 
include the effect on the albedo of clouds, atmospheric 
turbidity, relative humidity, and vegetation. The optical 
properties of these factors and their relationship to surface 
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temperature are less well known and cannot be easily para- 
meterized in a model as simple as the one at hand. 

The factor c in the outgoing infrared radiation 
term Rq , 

g 

(2d) c = 1 - m tanh (c^u ) , 

expresses empirically the "greenhouse effect", i.e., the 
screening by the atmosphere, in particular by the clouds 
in it, of infrared radiation from the earth, thus preventing 
part of it from reaching outer space. Notice that c decreases 
as u increases; this indicates that cloud formation, and 
hence the opacity of the atmosphere to low-frequency radiation, 
increases with increasing temperature. 

The function k((f>,u) in (Id) has the form 

C 4 — c s^ u 

(2e) k(4>,u) - k. [ ((l))+k 2 (4))g(u) , g(u) = — j e = f 1 (u) ; 

k^(<f>)u^ is sensible heat flux in the atmosphere and in the 
oceans, whereas k 2 (<f>) g (u) u^ is latent heat flux in the 
atmosphere. Here k n (4)) , k 9 (4>) are eddy dif fusivities , 
which parameterize convective transports; true conduction 
is known to be negligible in the atmosphere-ocean system 
on the planetary scale. In the original Sellers model, 
k(<f>,u) had the form 

k (<f) , u) = k g (<j),u) - v(u^) • (u+f (u) ) , 

where v is, for the purposes of modeling heat flux, mean 
meridional velocity. The theoretical shortcomings of the 
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additional terra v(u+f(u)) were pointed out by Robinson 
(1971) ; the practical difficulties in giving a good para- 
meterization are discussed by Sellers (1973). 

Numerical studies of Schneider and Gal-Chen (1973) 
indicate that results with the original Sellers model 
(denoted by them as (S) ) were very similar to those with 
the model adopted here (denoted in their work as (SV) ) , 
provided that the numerical values of k and k s were 
properly chosen (see the discussion on the determination 
of coefficients further on) . Furthermore, the recent 
work of Gal-Chen and Schneider (1975) shows that there are 
theoretical grounds oh which the formulation (SV) with zero 
meridional velocity is to be preferred. These considerations 
will be touched upon later, in a different connection. 

The constants c^ , 1 _< j £ 5, u^ , a , m , as well as 
the empirical functions C(4>), Q{<J>)/ b (<j>) , z(0)r k^((f)) and 
k 2 (4 > ) are made to fit currently measured values of temperature 
radiation, elevation, albedo and heat flux. The functions 
C(cj>) , Q(<{>) and z(0) are determined directly from measurements 
The function b(<j>) and the constant c^ in (2a) were 
computed by Sellers (1969) so as to fit existing albedo 
measurements in the northern and southern hemisphere. 

The constants m and c^ were also computed by Sellers, so 
as to fit empirical data on ; a is the Stef an-Boltzmann 
constant. The form of the function g(u) and the constants 
c 4 , c,_ appearing in (2e) are derived from certain physical 


- 9 - 



considerations having to do with the thermodynamics of wet 

air and from corresponding empirical data (see Handbook of 

Meteorology, 1945). The functions k 2 W are computed 

from measured data on sensible and latent heat flux, k. (4>)u, 

1 $ 

and k 2 ( <j>) g (u) u^ , respectively. These computations are 
based on the measured temperature distribution , denoted 
hereafter by 

u = u(<J>) , 

which will be called the data climate . Note that we use 
here the term "climate" only for convenience, instead of the 
lengthier "temperature distribution"? u({j>) is not necessarily 
a steady -state solution of the model? we return to this point 
in Section 4. 

The measured data are available at intervals of 10° 
latitude and are given in Table 1. Since the previously 
quoted authors used finite-difference formulations with a 
fixed 10° grid (except Faegre (1972) , who used a 5° grid) , 
these data were sufficient for their numerical work. In our 
numerical work, however, variable grid size was employed, 
and the 10° data were accordingly fitted by smooth functions. 

In fact, in order to have an additional check on the well 
posedness of the model (i.e., the continuous dependence of 
the solutions on the data, commonly referred to in the 
meteorological literature as sensitivity) , two forms of curve 
fitting were used: (i) by Bernstein polynomials, and (ii) by 

cubic splines. Results with the two forms of curve fitting 
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were indeed very similar (see Table 2). 

In this work we only consider symmetric solutions of (1) ; 
all data are symmetrized with respect to the equator, 

= tt/2. For such data the appropriate boundary conditions 
are 

(3a) u^CO) = 0 , (3b) u^ ( tt/2 ) = 0 ; 

in the symmetric case these are equivalent to u^(0)= u^(tt)= 0. 

We feel that the slight asymmetry between the northern 
and southern hemispheres could hardly have had a major 
influence on climatic change. Indeed, the circulations 
of the two hemispheres are practically separated from each 
other by the intertropical convergence zone, which acts with 
respect to our model as an insulating wall. The approximation 
involved in placing this wall at the equator is no worse, 
than other approximations in the model (see also Held and 
Suarez, 1974). A further reason for symmetrization will 
become apparent in the next section. 



3. Previous Results 


Budyko (1969) and Sellers (1969) used iterative numeri- 
cal techniques for constructing solutions of their time- 
independent models. They explored a range of values of the 
parameters appearing in the model, especially of the solar 
radiation Q,~~ and obtained one solution for each set of “ 
values. These solutions did- not depend smoothly on the 
parameter values; in several instances, small changes in the 
parameters led to large changes in the solution. For 
instance, an increase or decrease of a few percent in Q 
resulted in temperature changes leading to extensive melting, 
or significant expansion of the polar cap, respectively. 

Faegre (1972) obtained for a certain given set of values 
of the parameters five distinct solutions of his variant of the 
Sellers time-independent model. Two of these were highly 
asymmetric, and disappeared when c(<f>,u) in (1) was taken 
as constant? hence Faegre considered these solutions to be 
spurious and unphysical. It was the desire to eliminate a 
priori such solutions that suggested the choice of symmetric 
coefficients. Faegre 1 s formulation of the albedo is slightly 
different from that of Sellers, mainly in that the minus 
subscript in (2a) (i.e., the cutoff of a(<j),u) at u ) was 

missing. 

The three symmetric solutions of Faegre could be 
described as corresponding to the present climate, an 
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ice-age climate (about 15° C colder on the average than the 
previous one), and a completely ice~covered earth (at an 
average temperature of about 175 K) . This last climate 
was also the one obtained by Sellers when decreasing the 
solar radiation by more than 4%. 

These results raised the question of the transitivity 
of the earth's climate, as formulated by Lorenz (1968, 1970). 
In Lorenz's terminology, a time-dependent system of equations 
is transitive if its solutions have a unique equilibrium 
statistic, that is, if all solutions, independently of 
initial conditions, have the same infinite time average; 
otherwise the system is intransitive . Lorenz pointed to tne 
existence of certain transitive systems which possess a 
property called by him almost intransitivity , i.e., that 
of having at least some solutions whose averages over long, 
but finite , time intervals are different — these solutions 
then would alternate in time between the different averages. 
He raised the possibility that the atmospheric system is 
almost intransitive? in other words, that the known climate 
changes in the past were not necessarily caused by changes 
in external conditions (like solar radiation) , but rather 
were an effect of the normal evolution of the system. 

Schneider and Gal-Chen (1973) investigated the question 
of transitivity for energy-balance climate models by 
formulating time-dependent versions of the Budyko (B) , 

Sellers (S and SV) and Faegre (F) models. They solved 
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numerically the initial-value problem governing these time 
dependent models for a large range of initial conditions. 

The models were found to be intransitive, rather than almost 
intransitive: every solution tended as t -*■ 00 to one of 

two (or more, in the case of the (F) model) equilibrium 
solutions; viz., the equilibrium statistic of the system 
governing each model was not unique, and no spontaneous 
transition from one equilibrium to another was possible. 

The two equilibrium solutions obtained for all models 
corresponded to the present climate and to the previously 
mentioned completely ice-covered earth. 

The equilibrium solutions, at least for the (S) and 
CSV) models, proved stable under rather large perturbations 
in both initial conditions and parameters. That is, solu- 
tions which differed in their initial conditions from one 
of the limiting ’’equilibria" by as much as + 18 K tended 
as t -*■ 00 to the "equilibrium" near which they started. 

Also changes of +1.5% in the solar radiation led to limit- 
ing equilibria which differed by only a few degrees from 
the unperturbed ones. However, changes of more than - 18 K 
in initial conditions or - 1.5% in solar radiation led from 
the. "present climate" to the "ice-covered earth". The latter 
seemed to be the most stable climate in all investigations 
mentioned above . 

Schneider and Gal-Chen did not obtain a limiting steady- 
state solution which would correspond to a true ice age as 
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recorded in the planet's history. Their results with the 
time-dependent version of the Faegre model (F) were rather 
different from those with the two versions of the Sellers 
model (S and Sv) , especially with regard to the stability 
of steady states under perturbations. 

Contrary to the results of Schneider and Gal-Chen, 

Dwyer and Petersen (1973), with a time "dependent model 
essentially identical to Schneider and Gal-Chen's (S) , 
obtained solely one type of limiting steady state/ that 
corresponding to the "present climate". They used only the 
data climate u = u(<J>) as initial conditions, but varied 
the solar radiation Q. The actual values of the limiting 
equilibrium depended of course on the values of Q used, 
but slight changes in Q yielded only a difference of a few 
degrees between the average temperature of the equilibrium 
and that of the data climate; no "deep freeze" equilibrium 
was obtained. 

These results seemed to be interesting enough in order 
to warrant further study of energy-balance climate models. 

As indicated in the previous section we chose to investigate 
symmetric solutions of the (SV) model of Schneider and 
Gal-Chen, and of some variations thereof, including the (F) 
model. We hope that . this study will add as much light and 
as little heat as possible to the climate question (Jackson, 
1962). 
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. Steady-State Solutions 


We turn now to the mathematical theory of equation (1) . 
Introducing the new space variable x = 2<J>/7r, we obtain the 
initial-and-boundary value problem (IBVP) 



(5a) u (0#t) = u (1^ t) = 0 , (5b) u(x,0) = u(x) , 

A A 


where (4) is a nonlinear parabolic partial differential 
equation (PDE) . Here g(u) is given by (2e) , u - 1 
(its significance will appear later, in Section 6) and 

c 1 = 0.009, c 2 = 0.0065 deg irf 1 , c 3 ~ 1.9 x 10~ 15 deg" 6 , 
(6) c^ = 6. 105x0. 75xexp(19. 6) xio^dyn deg cm 2 , c^= 5350 deg, 

a = 1.356x10 ^cal cm 2 sec ‘'“deg ^ , m = 0.5, 283. 16 deg. 

Mesh data at 10° latitude for C(x), k^(x), k^fx), Q(x), b(x), 
z(x) are given in Table 1. The units of the constants and 
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mesh data are chosen such that the common units of all terms 
in (4) are cal cm ^sec \ 

The first step of the investigation is to find steady- 
state solutions of (4, 5a) in the range of physical interest, 
and to study their dependence on changes in the model. We 
consider therefore the steady-state equation obtained from (4) 
by setting = 0. After some rearrangement we get the 

following two-point boundary-value problem (BVP) for the 
system of ordinary differential equations (ODE) : 


(7a) 


u = v 

X 


(7b) 


v = - 

X 


rIL'i 2 F (x,u) 
k 2 } k (x, u) 


2~(cot 


TTX> 

2 J 


v - 


kj (x)+k 2 (x) g(u) 


k (x,u) 


v 


k 2 (x) g' (u) 2 

k (x, u) V 


0 < x < 1 


( 8 ) 


v (0) « v(l) - 0 , 


where 


(7c) 

(7d) 


k (x,u) 
F (x, u) 


= k^x) 

= yQ (x) 


+ k 2 (x)g(u) , 

|l - b (x) + c 1 [u m + (u-c 2 z (x) -u^) 

4 6 

- au [1-m tanh (c^u ) ] . 



We use shooting (Isaacson and Keller, 1966, Keller, 1968) 
as the numerical procedure to solve (7, 8) : equations (7a,b) 
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are solved with initial conditions 


(9a) v (0) = 0 , (9b) u(0) = u Q , 

for different values of u Q ; denote by u(x?u Q ), v(x;u Q ) 
the solution of the initial-value problem (IVP) (7,9) in 

order to emphasize its dependence on the parameter u Q . 

An iterative scheme is then used to obtain those values 
of u Q which satisfy 

(10) v ( 1 ; u Q ) “ 0 . 

For these values of u Q the solution u(x;u Q ), v(x;u Q ) 
of the IVP (7,9) is also a solution of the BVP (7,8). 

To obtain numerical solutions of prescribed accuracy to 
the BVP (7,8) one has therefore (Keller, 1968) to achieve 
the desired accuracy both in 
(i) solving the IVP (7,9), and in 

(11) solving iteratively the nonlinear (finite) equation (10) 
In solving numerically the IVP (7,9), we encounter the 

difficulty that (7b) is singular at the origin, since 
cot ( ttx/2 ) + “ as x 0. This singularity arises from the 
mathematical form of the diffusion term D of (1) in spherical 
coordinates . It is easily shown , though, that (analytical) 
solutions of the IVP exist and are unique at least "in the 
small", i.e., for | u Q -U 0 | < e , 0 < x < 6, arbitrary U Q 
and some 6, e > 0. The numerical difficulty of the initial 
point being singular can be circumvented by using a variable- 
step method. 
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It also turns out that in order to achieve prescribed 
over-all accuracy of the numerical solution of (7,9) with 
a minimum of computational effort it is convenient to use a 
variable-step variable-order multistep scheme. The ODE 
solver we computed with was developed at the Lawrence Liver- 
more Laboratory and is documented by Hindmarsh (1972,1974). 
Both the Gear and Adams methods included in the package were 

used and gave results identical to within the prescribed 

-7 

relative accuracy, which was 10 ; however, the Adams method 

was faster and was therefore used in most integrations. 

The number of steps needed per solution (per value of u^) 

for this accuracy was of the order of 100. It is of interest 

to point out that near the singularity at x = 0, and near 

the point where the jump discontinuity in the derivative of 

F(x,u(x?Uq)) occurred, the order of accuracy chosen by the 

scheme was one and the step size was very small, so that- a 

proportionately larger amount of computation was done in the 

neighborhood of these two points. 

Thus every solution of (7,9) for given u^ puts at 

our disposal a value of the function v(1;Uq), accurate to 
-4 * 

10 . To find the zeros of v(l;u Q ), we used the method of 

false position (Isaacson and Keller, 1966). The criterion 

— O 

for Uq to be a root of (10) was |v(1,*Uq) | < 10 . 

** ' 3 

In the units chosen, typical values of u are 0(10 ) and 

2 

of v are 0(10). 
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The range of physical interest in which we searched for 
solutions of (10) was 100K £ u Q £ 300K (see also Faegre, 
1972). Within this range, three solutions of the BVP (7,8) 
were found which are denoted by u^x), u 2 (x), u^(x). 

They correspond to 


(0) = 247.74 K , u 2 (0) = 223.97K u” (0)" "= 168.94 K. 


The curve u (1) vs. u(0) , the zeros of which correspond 
X 

to the solutions u 3 ' u 2 ' u i ' is 9 iven in Figure 2. The 
individual solutions u-^(x), u 2 (x) , u^(x) are plotted in 
Figure 3. 

It is customary and useful to characterize a climate 
of the model, i.e., one of the solutions above, by its 
average temperature, rather than by the temperature at the 
pole. Therefore we introduce for functions <f> the 
averaging operator £ by 


£ cj> 


-f 


( x) sin 



jx, , 

sxn (— ) dx . 


In particular , it is known (Isaacson and Kel ler , 1966) that 
for functions <J> symmetric about the equator, <f>(2-x) = (x) , 

which also satisfy <J> (0) =0, and hence can be extended so 
as to have period 2, the trapezoidal-quadrature approximation 
of £ is very accurate. Denote this approximation by £^ , 
where A = 1/9 (corresponding to a 10° mesh) ? then we have 

£ a u x = 287.76 K, & A u 2 - 267.44 K, & A u 3 = 175.43 K . 
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Clearly u^ corresponds to an interglacial, to a glacial, 
and u^ to the ice- covered earth. 

Note that for the data climate u = u(x) 

u ( 0 ) = 247.36 K, = 287.20 K , 

which is indeed very close to the interglacial solution of 
the model, (see also Figure 3) . Though present observa- 

tions were used to obtain constants and empirical functions 
in (4) , no explicit term was added to ensure that (4) hold 
with u fc = 0; also the diffusion term, D, is of the same 
order of magnitude as the radiation terms, and Rq , so 
that this is not the result of a simple radiation balance. 

The same is true of the work of all the previous authors 

* 

discussed in Section 3 ; these authors, however, assumed 
at least implicitly that the data climate should be a 
steady-state solution of the model, or approximately so. 
Since this result is not actually built into the model, as 
far as we can see, we think it is rather remarkable that 
this class of models yields a steady-state solution, u-^(x), 
very close to the data climate, u(x) . Henceforth we shall 
refer to also as the "present climate" of the model. 

* In fact, Schneider and Gal-Chen (1973) did use an extra 
"fudge coefficient", 0.97 < c Q (x) £ 1.03, in R Q in order 
to achieve better agreement between the interglacial, or 
"present" climate of their models and the data climate; 
they obtained &u^ = 287.06 K vs. £u = 287.30 K for the 
( S ) mode 1 . 


- 21 - 



The mesh data, from which constants and empirical 
functions for the model were computed, are known to be in 
error by possibly as much as 100% (Schneider and Gal-Chen, 

19 73); also some of the parameterizations used in formulating 
the model are questionable. Therefore we made a number of 
com putati ons in order to obtain information on the dependence 
of the physically significant solutions, u-^u^u^, on varia- 
tions in the model. The salient features of these computa- 
tions, summarized in Table 2, are: 

(1) The dependence of all these solutions on the 
coefficients seems to be very smooth. 

(2) The bounds on the albedo, 0.25 £ a(x,u) £ 0.85, are 
essential for the existence of three steady-state solutions 
in the physical range, 100 K £ u(0) £ 300 K: u 3 disappears 
when the bound a £0.85 is not enforced, and u^ disappears 
when the bound 0.25 £ a is not enforced. 

(3) The terrestrial radiation term Rq , given by (lc), 

and its nonlinearity are essential for the existence of 
physically significant solutions: if the term is set equal 

to zero all solutions, u^,u 2 ,u^ , disappear; when it is 
replaced by its average, 

~ ~ 4 

r q = (x,u(x) ) au (x) = const. , 

u 2 only obtains; when Rq is replaced by a linear approximation, 

R q = au 3 {u+4(u-u) }£ a c(x,u(x) ) , u = £^u(x) , 

then u^ only obtains. 
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(4) Our results with the model using the Faegre formula- 
tion of the albedo are very similar to those using the Sellers 
formulation; it is hard to explain the disagreement in this 
respect between our results and those of Schneider and 
Gal-Chen (1973) , who, as mentioned in Section 3, obtained very 
different results when using the two albedo formulations. 

(5) The values of the empirical functions in the model 
when fitting mesh data by (i) Bernstein polynomial approxima- 
tion and (ii) cubic spline interpolation are pointwise rather 
different for some of the functions (see Figure 1) . The 
solutions of (7,8) when using these different fits are, however, 
very close. This also supports the assertion in paragraph (1) 
above, and shows that the uncertainty in the mesh data does 

not affect in an important way the conclusions of the investi- 
gation. 

We explored another variant of the model, in which the 
eddy diffusivity k(x,u) given in (7c), corresponding to the 
(SV) model of Schneider and Gal-Chen, was replaced by k(x,u), 
where u is the data climate; hence also g’(u)v in (7b) is 
replaced by g' (u)u' (x) . This variant, in accordance with the 
terminology of Gal-Chen and Schneider (1975) , would correspond 
to an (SVC) model'. We denote this modification of (7.) by (7'); 
the solutions of (7',8) are 

u L (0) = 247.55 K, u 2 (0) = 227.76 K, u 3 <0) = 169.44 K, 

& A u 1 = 287.70 K, & A u 2 = 268.60 K, & A u 3 = 175.44 k. 
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These are indeed very close to the solutions of (7,8), 
especially for u^ and (see also Figure 3) . 

Gal-Chen and Schneider (1975) investigated the effect of 
variations in the solar radiation yQ on the equator-to-pole 
temperature gradient; they argue that this dependence should be 
monotone, in fact monotone increasing- In the light of this 
argument and of their results, model (SV) is superior to (S) , 
as already mentioned in Section 2; moreover, model (SVC) is 
superior to (SV) . Since (SVC) is also more convenient to use 
when investigating the stability of the steady-state solutions, 
we shall work with it in the sequel; the corresponding form 
of (4) we denote by (4' ) , the corresponding form of (1) by (1*), 
and the corresponding form of (7) by (7'). 

For this model, additional computations were made outside 
the range of physical interest. One more steady-state solution 
was found; we denote it by u 4 (x) , and it satisfies 

u 4 (0) = - 185.99 K, & A u 4 = - 175.40 K . 

Most probably there are no other solutions of the BVP (7',8) 
at all. Indeed, the solutions of the IVP (7* ,9) become 
unbounded as Uq moves towards the ends of the range explored, 
viz., -1330 K £ u Q < 300 K; i.e., u(x;u Q ) ■* +°° as u Q approaches 
the ends of the interval above (see Figure 2). 

We would like also to point to the fact that some of 
the modifications of the model which have three physical 
steady states yielded numerical values of the latter consider- 
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ably higher than u-^u^u^ (see Table 2); however, it was the 
qualitative feature of the number of solutions and their 
approximate position with respect to each other that was of 
interest in our investigation: the average temperature of 

a given solution of any fixed model could be changed to 
practically coincide with that of the solution of (4) to 
which it corresponds by adjusting the values of the coeffi- 
cients. 
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5 . Stability of the Steady-State Solutions 

In the previous section we have shown that equation (4') 
with the boundary conditions (5a) has three steady-state 
solutions of physical interest, u = u^(x), 1 £ j £ 3* In 
this section we shall study the stability of these steady 
states. 

Let us concentrate on any one of the steady states above, 

u = u.{x), where i is fixed. Stability of u. means that small 
3 3 

perturbations of u^ die out with time. More precisely, 

o 

is stable if, when taking any nearby state U, i.e., one 

which differs little from u. , 

3 

(11) U = Uj (x) + ev(x) , 

o 

say, where v is arbitrary and e > 0 is not too large, then 
the solution U(x,t;e) of ( 4 * — 5 ) with initial condition 

o 

U(x,0;e) = U(x;e) tends to u^ itself as t + » . 

Let (4') be written symbolically as 

(12) u t = N (u) , 

where we divided through equation ( 4 1 ) by C (x) , 0 < C m <_ C(x) 

£ C M < «> , and N(u) is the corresponding right-hand side. Take a 
perturbed u, u = U(x,t,e), which is assumed to satisfy 
the boundary conditions (5a) and the initial condition (11) . 

If such a solution of (12) exists for all e sufficiently small, 
0 < e 1 £ q • then the equations 
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( 12 ') 


(x,t;e) -U t (x,t ; 0) 


N ( U ( x , t ? e ) ) -N (U (x, t;0) ) 


0 < e < £ Q , 


also hold, with U(.x,t?0) - (x) . Letting £ + 0 we 


obtain the linearized equation 


(13a) 


v t = - L.V , 
_ 3 


where we define v(x,t) = U(x,t;e)| £ __ and 


(13b) L. = 


f T N(u(x,t;e ))| £=0 = 


- t— N (u) [ , 1< j <3 

3u u=u . — J — 


At this point, for the sake of simplicity, we shall 
drop the subscript j, and u will thus stand for some u^ , 
L for the corresponding L ^ . 

Equation (13a) is linear in v and it has a unique solu- 
tion v satisfying the boundary conditions 


(13c) 


v ( 0 , t) = v ( 1 , t ) = 0 , 

A X 


and arbitrary initial conditions 


(13d) 


v(x, 0) = v ( x) 


The solution may be found by the usual method of separation 
of variables, or expansion in eigenfunctions (normal modes) 
Consider 


v(x,t) = e ^ w(x) ? 
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this v will be a solution of (13) with v(x,0) = w(x) iff 
(if and only if) A is an eigenvalue and w is an eigenfunction 
of L, i.e., iff the pair (A,w) satisfies the homogeneous 
problem 

(14a) Lw = Aw 

(14b) w (0) = w (1) = 0 . 

X X 

We shall show that the operator L has sufficiently many 
eigenfunctions, and that therefore any solution v of (13) 
can be written as a series 

” .iOO*. (V) 

(15) v ( x , t ) = l cLe w { } (x) 

k=0 K 

(k) (k) 

where (A ,w' ') are all the solutions of (14). It is 

clear from the derivation ( 12 1 ) of (13) and from (15) that 

for the steady state u to be stable it is necessary that 

(k) 

every eigenvalue A of L have positive real part: 

(16) Re X (k) > 0 . 

This condition defines the so-called linear stability of u. 

It has been shown rigorously in a few cases and it is believed 
in many others that (16) is also sufficient for the stability 
of u as we defined it at the beginning of this section, i.e., 
that linear stability implies nonlinear stability. In the 
next section we shall give an argument which will make it at 
least plausible that this is the case also for the problem 
at hand. 
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We turn now to the analysis of the eigenvalues of the 
linear second~order ordinary differential operator L, which 
we can write as 


(17a) 

LW as - 

’ r(x) ^ p(x)w x^x + ‘I (x)w • 

Here p 

,q,r are determined by (13b) and by (4' 

(17b) 

r(x) = C(x) sin (-^) , 

(17c) 

p(x) 

= (^ ) sin ( 2 ) * k (x, u) , 

( 17d) 

q (x) = 

{Q(x)(c 1 ) c , - d (x, u) }/C (x) , 

with 


f c^ if .0.25 < a(x,u) _< 0.85 

( 17e) 

(c l>c' = ' 

and u (x) -c^ z (x) -u < 0 
2 m — 



^ 0 otherwise , 

and 



( 17 f ) 

d (x, u) = ~ c(x,u)au 4 


= au^|4[l-m tanh(c 3 u 6 )]- 6mc 3 u^ [ 1-tanh 2 (c^u^ ) ] j- 

Clearly L thus defined is formally self-adjoint (Courant 
and Hilbert, 1953, Birkhoff and Rota, 1969) under the 
prescribed boundary conditions with respect to the inner 
product 

1 

(18) (w 1 ,w 2 ) = | r (x) w x (x) w 2 (x) dx . 

0 

Moreover, r,p are nonnegative and twice continuously differ- 
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entiable on the interval I = [0,1] r and q is piecewise contin- 
uous on I. However, because of the singularity at x = 0 due to 
the fact that r(0) = p(0) - 0, the usual Sturm-Liouville theory 
of self-adjoint operators (Courant and Hilbert, 1953, Birkhoff 
and Rota, 1969) does not apply to L; this difficulty though 
-can-be._ove.rcome_and_ the - theory can be... extended.. . The crucial - 
element in this extension is the observation that, because 
of the boundedness of q, L is bounded from below in the sense 
that, for any w satisfying (14b) for which (w,w) < 00 , the 
inequality 

(19a) <w,w> £ K(w,w) , 

holds with some fixed constant K, K > min q(x) >-», independent 

0 <x<l 

of w; here <*r,w> is the Dirichlet integral corresponding to L, 

1 

f 2 2 

(19b) <w,w> E (Lw, w) = (pw x . + rqw ) dx . 

0 

This result, together with an analysis of the nature of 
the singularity at x * 0, are sufficient to show that indeed 
L has a complete system of eigenfunctions, orthonormal with 
respect to the inner product (18) (Birkhoff and Rota, 1969) , 
and that the eigenvalues of L are real and can be arranged 
in an ascending sequence 

-«■ < A (1) < A (2) < A (3) < ... 

with A^k^ -v oo as k *> °°, and K = A ^ in (19) . Hence 
any solution v of (13) can be written as a series (15) . 
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Therefore the stability question for the steady state u 
reduces to determining whether 

(16') A (1> > 0 , 

in which case u is stable , or whether the opposite holds, 
in which case u is unstable. 
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5a. Stability Criterion 


In this subsection we shall show that it is possible to 
determine the sign of A^ without actually computing A^. 

We start with the variational characterization of the 
lowest eigenvalue (Courant and Hilbert, 1953), A^, 

(2 0) A ^ = min = min <v,v> , 

<v , v> <°° lv ' v; (v,v)=l 
(v, v)^0 

i.e., A is the minimum of the Rayleigh quotient R(v) 
corresponding to L, 

(21) R ( v) = <v,v>/(v,v) , 

and the minimum is assumed for v = w^. Indeed, (14a) is 
the Euler equation for (20), and v (1) = 0 is the "natural" 
boundary condition for (20) in the sense of the calculus of 
variations (Courant and Hilbert, 1953) ; also v (0) = 0 is, 

X 

according to the theory of singular operators with the 
properties of L, the only boundary condition at x = 0 which 
ensures that, for solutions v of (14) , <v,v> as well as (v,v) 
is finite. In particular, we conclude that the minimizing 
function, v = w^ , is at least as smooth as the coefficients 
of L, viz., it must have a piecewise continuous second 
derivative. 

With these preliminaries we are able to prove the 
following known 
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Lemma . The first eigenfunction of L, w^, is strictly 
pos i ti ve , 

w (x) > 0 , 0 1 x 1 1 • 

Proof : From the definition (21) of R(v) it is clear that 

R(|w U) |) = R{w (1) ) , 

where |y| denotes the absolute value of y. If w^ had a 
zero at some interior point x Q , 0 < x Q < 1, and (x Q )^0, 

then the first derivative of |w^ | would have a jump at 
x - Xq , which contradicts the smoothness of |w^ | as a 
solution of the variational problem. If, on the other hand, 
w^ 1 * (x 0 ) =0, 0 < x 0 < 1, or w (1) (0) = 0, or w (1) (1) = 0, 
then, by the uniqueness theory of linear ODE (Birkhoff and 
Rota, 1969), it would follow that w^ = 0; this contradicts 
the fact that is a nontrivial solution of (14), 

completing our proof. 

Now we are ready to state our stability criterion, which 
is part of the conclusion of the 

Theorem . Let L, A = A^^, w = w^ be as above. Suppose 
also that there exists a function v, as smooth as w, 
nonnegative, v 0 on I , and satisfying 

(22) Lv _> 0 , v (0) = 0 , v (0 ) = 1 . 

X 

Then v (1) 0 implies * >_ 0 . Moreover, 

(i) A > 0 if either Lv % 0 or v (1) > 0, and 

X 

(ii) if Lv = 0 holds, then A>0, A-0, or A<0, according to 

whether v (1) > 0, v (1) =0, or v (1) < 0. 
xx x 



Proof ; The required results are easily read off from the 
following sequence of eq -ities obtained from the definitions 
and by integration by parts; 


X 

X j rvw - 
0 


A(w,v) = (Lw,v) = I -(pw ) v + rqwv 

J X X 

0 

1 .. 1 


-pw v 

X 


+ pw v + rqwv 

J X X 

0 o 

1 l 


pwv 


X 


+ J “ (P v x ) x w + rqvw 
0 0 


= (pwv ) (1) + (Lv,w) , 

A 


where we used w (0) = w (1) = 0 and v (0) ~ 0 . 

X X X 

Note, This is essentially a comparison theorem, of the 
type familiar from Sturm-Liouville theory (Birkhoff and Rota, 
1969). 


The result under (ii) above is the stability criterion 
which we shall use. For completeness we give here also the 
following slight generalization as a 


Corollary . Let L, X, w be as in the Theorem. Suppose v 
satisfies the hypotheses of the theorem, except for that of 
being nonnegative on I. Instead, let x^ , 0 < x^ £ 1, be 
the first zero of v, 

v(x) > 0 on 0 £ x < x Q . 

Then X < 0. 
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Proof : The result follows from the same integration by 

parts carried out in the proof of the theorem, except that 
now the upper limit of integration has to be x^ rather than 
1, and we use v(Xq) = 0 rather than w^(l) = 0 in passing 
from the second to the third line. Moreover, since v(x) > 0 
for x < Xq , and v is continuously differentiable, v^(Xq)£ 0. 
Furthermore, it cannot be that both Lv = 0 and v^(Xg) = 0, 
since then we would have, by uniqueness, v = 0, which 
contradicts v(0) = 1. This completes the proof. 
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5b . Stability Results 


In this subsection we shall apply the stability criterion 

obtained in Subsection 5a to the steady -state solutions u. , 

3 

1 £ j £ 3 * 

For this purpose we construct functions by solving 
(22) with the equality sign, and with L = given by (13b) 

and (17). The results, obtained with a relative numerical 

-7 

accuracy of 10 , are that 

(a) Vj , j - 1,2,3, is positive, v j £ V , > 0, >_ 0 . 5 ; 

(b) Vj(l) > 0 for j = 1,3, whereas ^ v^ (1) < 0 for j=2. 
The actual computed values are 


v, (1) = 2. 39970, 
dx 1 


= 


-3.24312, ^v 3 (1)=4. 31025 


These values of the derivatives at x = 1, together with the 
values of , are sufficiently bounded away from zero in order 
to conclude that the conditions of the theorem are satisfied 
beyond the doubt of numerical uncertainty. 

Hence the solutions representing the present climate and 
the ice-covered earth are stable, whereas the so-called ice 
age of the model is unstable. This result agrees with and 
throws additional light on the results of previous authors, in 
particular the time -dependent integrations of Schneider and 
Gal-Chen (1973) . 
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Actually the stability computations were carried out 
also for the solutions of 


(4") 


N ' { u j ) - 0 , 


1 < j < 3, 


with 


(23) N' (u) = (r) 2 


K n o 

0 d . ,TT X, 

tt / sin(77x/2) 3x Sin 2 u x 


+ Q(x)|l - B 0 + c^u 


c “ C 6 au 


where the subscript 


{}. 


is defined in (2c) and 


K Q = S A k(x f u(x)) * 2.2x10 5 , B Q = (x) * 2.85881, 

c 6 = £ a c ( x 'U(x)) = 0.61 ; 


the corresponding linearization of N 1 is 

(13') l'. = - N (U.) = (-) 2 |- 

j u 3 IT sm(Trx/2) 3x 


sin (— } H 


where 


+ Q(x)(c 1 ) c „ - ^ c g au j , 


(C 1 J c" 


c. if 0.25 < 1-B-+C.U. < 0.85, 

1 — 0 1 j — 

0 otherwise 


This is, according to the experiments carried out in Section 4 
(see Table 2), the simplest form of (4) which still yields 
approximately the same steady -state solutions in the physical 
range of interest. The results are the same, i.e., u^ and u^ 


- 37 - 



are stable, and u 2 is unstable. It seems therefore quite 
plausible to conclude also that for all models, lying in some 
sense between (4) and (4"), the steady states corresponding 
roughly to u^u^u^ have the same stability properties as the 
latter. Combining this remark with the one made at the end of 
Section. A ,__.i_t seem s that w e have a . result. . about, the . s tabi 1 i ty 


of the steady states for a certain type of energy-balance model, 
rather than just for one specific model of this type. It seems 
desirable to define precisely the type of model having these 
properties, and we intend to try to do so in further work. 

Having thus determined the stability properties of the 
steady states of (4’),we want to obtain some additional infor- 
mation by actually computing the lowest eigenvalues a!^ r an< ^ 
corresponding eigenfunctions * 1 1 j 3. The eigen- 

values are 


A^ 1 * = 3.95390x10 9 


( 1 ) 

2 


-9 


AA~' = -5.87662x10 ", X;*'- 5.13587x10 , 


(U_ 

3 


-9 


( 1 ) 


and the eigenfunctions wj*' , 1 £ j £ ^ are plotted in Figure 4. 

The method for computing was again shooting, 

this time with respect to the parameter A. That is, equation (14a) 
with L = Lj , was solved for w(x;A) with initial conditions 
w (0; A) = 1 , w (0;A) = 0 , 

X 

and with different values of the "shooting parameter" A. The zeros 

of the function v x (l;A) were then computed by regula falsi to 

-4 . 

within an accuracy of 10 . The over-all relative error m 

-7 

computing solutions of the initial-value problem was 10 , 

as before. 
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Notice from (15) that 


T = 1/| Xj 1 * | 

is a characteristic decay or relaxation time for the solution 
U(x,t;e) of (12) to u = u^ (x) . This time is of the order 
of 10 years for all j , 

= 8.0 years , t 2 = 5.4 years,, = 6.2 years. 

In this context it is remarkable that Schneider and Gal-Chen 

(1973) state that, in one of their integrations, the solution 

o _ 

of (4) with initial data u(x) = u(x) and y - 0.984, after 

dropping rapidly by about 12 K in average temperature, was 

nearly constant f or . about 50 years of simulated time, and 

then finally dropped to a steady state close to our u^(x) . 

From our results it becomes clear that the solution mentioned 

above hovered for a time of the order of t 2 around u 2 r but 

could not persist there indefinitely because of the instability 

* 

of u 2 , and finally attained u 3 , which was stable. 

It also follows from (13) and (14) that multiplying C (x) 
by a constant k > 0 will result in Tj, 1 < j < 3, being 
multiplied by k. A similar statement holds also for the 
nonlinear equation (4), since such a constant k just 
corresponds to a different scaling of the time t (see also 
Schneider and Gal-Chen, 1973). Experiments in determining 

* We shall see in the next section that x 2 increases with 
decreasing y. 
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characteristic response times for solutions of (4') were 
done by Dwyer and Petersen (1973) , who used two heat capaci- 
ties C(x) , both larger than that used by Schneider and 
Gal-Chen (1973) and by us. It seems, however, that upon 
decreasing y in (4) to y = 0.98, as they always took 

o 

u(x) = u(x) in (5b) , the average temperature of the^ solution 
u(x,t) dropped rapidly at first, and only slowly thereafter, 
as indicated by Figure 2 in their article. Apparently this 
slow decrease, which shows that their solution u(x,t) of (4 1 ) 
was approaching a steady state close to our U 2 (x), was 
interpreted by them as proving the nonexistence of a steady 
state u^(x) , which had been obtained in the work of Budyko, 
Sellers and Faegre when decreasing y by a similar amount. 

The influence of changes in y on the steady-state 
solutions Uj , 1 < j < 3, of (4') will be investigated in 
the following section. At this point we only want to mention 
that, whereas a multiplicative factor k in C(x) affects 
the magnitude of and hence of T j' our theorem shows 

that it does not affect the actual stability of u^ , i.e., 
the sign of • 
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6 . Perturbed Steady-State Solutions and Bifurcation 


It is clear that steady-state solutions of (4') could 
not spontaneously evolve into each other. More precisely, 
the solution of (4'-5a) with initial condition u(x,0) = u^ (x) , 
j = 1,2,3, is u ( x , t ) = Uj (x) . Moreover, it stands to 

o 

reason that, for any physical initial condition, u(x) _> 100 K, 
we would have 

lim u(x,t) = u. (x) , 

t~K» 3 

with j = 1 or j = 3, since u 2 (x) is (linearly) unstable, 
whereas u^(x), u^ (x) are (linearly) stable (and u^ (x) <-170K 

< 0) . 

Thus, to explain physical ice ages, one has to consider 
perturbations in the parameters appearing in equation (4‘). 
Such perturbations would presumably be caused by physical 
mechanisms not included in the model. The most likely 
candidate for such a role is ]i , which up to now was taken 
to be unity. Indeed, many ice-age theories rely heavily on 
a change, however small, in the amount of solar radiation 
reaching the lower layers of the atmosphere (SMIC, 1971, 
Beckfnsale, 1965) . 

Some attribute the assumed decreases in solar radiation 
to changes in the parameters of the motions of our planet 
(Milankovitch, 1930) , others to airborne volcanic dust 
due to an increase in volcanic activity (Fuchs & Patterson, 
1947), and so on. There has also been concern about a 


- 41 - 



possible climatic catastrophe being imminent because of the 
increase in the quantity of industrial pollutants in the 
atmosphere (Rasool and Schneider, 1971) . 

To investigate the effect of such changes in the model 
at hand, the curve u (1) vs. u(0) was recomputed for 
di f ferent_ values._of __y, in particular with a view to obtain- 
ing u^(x;y), u 2 (x? y) . One important result is that these 
two solutions coallesce for 

y c « 0.98152 

and disappear entirely for y < y . The bifurcati ng solution 

c 

u = u (x) = u.. (x;y ) = u 9 (x?y ) was computed with over-all 
C X c Z C 

— 9 * — 7 

relative accuracy 10 and |u (1)| < 5xio . The 

(u(0) ,u (1) ) -curve corresponding to y = y is given in Figure 5 
X c 

notice that it is very flat near u(0) = u (0) , which makes 

c 

it difficult to compute U Q ( X ) accurately. 

Further computations were carried out for y = 0.982, 

0.985, 1.01, 1.02, 1.03, 1.04 and 1.05. The results of these 

computations are summarized in Table 3 and plotted in Figure 1 6. 

It is quite interesting that, whereas for y _> 1.0 the 

dependence of u^ (0) , £^u_. (x) , j = 1,2, on y is almost linear, 

this dependence is definitely quadratic near y = y ; the 

c 

latter is in good agreement with the general theory of 
bifurcation for nonlinear parabolic problems (Hoppensteadt 
and Gordon, 1975). According to the theory, there exists 
in principle the possibility that, instead of disappearing 
at y = y c , the two solution branches u^(x;y), u 2 (x?y) could 
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merge into one periodic solution for y < y . 

This possibility was not borne out, however, by the time- 
dependent computations of Dwyer and Petersen (1973), and of 
Gal-Chen and of Schneider (1973, 1975) ; we did not investi- 
gate it further. 

Concerning the problem of the pole- to- equator tempera- 
ture gradient, 

AUj = u j(D ” u^(0) ' 1 i j £ 3, 

as discussed by Stone (1973) and Gal-Chen and Schneider (1975) , 
the curve Au_. = Au^ (y) , j = 1,2, is particularly interesting. 
We notice that 

(a) the values of Au-^ lie below those for Au 2 ? 

(b) the values of Au^ are monotonically decreasing with y; 

(c) the values of Au 2 have a maximum for y somewhere between 

y = 0.985 and y = 0.99; 

(d) the dependence of both Au^ and Au 2 on y, but especially that 
of Au, , is very steep near y « y . ■ 

-L c 

Being aware of the limitations of the model, as pointed out 

in Section 2, ,we do not want to make extensive comments 

concerning these results, but only note them for comparison 

with the results of other models and for further study. 

We also studied the stability of the solutions u = u c (x) 

and u = u^.(x;y ). Repeating the construction indicated at 
j c 

the beginning of Subsection 5b for v = v c and v = v^ , where 
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L . 
3 


i 


j c f 3 , 



we obtained the following results: 

(a) Vj / j = c,3, is positive, v j — V j > V j 0.75 ; 

(b) t- v (1) = - 7.56X10" 3 , — v,(l) = 4.23362. 

Clearly u 3 (x?u c ) is still stable, in fact (d/dx) v^ ( 1; jj ) 

=4.23 is very close to (d/dx) v^ (1 ; 1. 0 ) = 4.31, showing 

the extremely smooth dependence of u^(x;ia) on u- 

The negativity of (d/dx) v (1) would seem to point to 

outright instability of u^ (x) , but in fact its small 

numerical value indicates that, within the accuracy of the 

computations, it is actually zero, i.e., u (x) is neutrally 

c 

stable. Indeed, the mathematical theory of nonlinear 
problems (Nirenberg, 1974) shows that for a bifurcation to 
exist, as in the case at hand, the linearization 


= - ^N( U ;uJ 


u=u 


of the spatial, time-independent, operator N at the bifurca- 
tion point (u ,u ) has to have a zero eigenvalue. This 

G G 

follows from the infinite-dimensional generalization of the 
one -dimensional fact that a function can be inverted in a 
neighborhood of a point at which it has a nonzero derivative / i »s., 
that it has a unique branch near such a point. 

We also computed the lowest eigenvalues 


x! 1 ’ and 


corresponding eigenfunction w 


( 1 ) 


of L j (]i c ) for j = c , 3 , 
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by the shooting method described in Subsection 5b. In this 
computation, the results on v ^ mentioned before were valuable 
in making a first guess for the shooting parameter A. The 
computations yielded 


A (1) = - 1.287x10 11 , A* 1 * = 5.07265x10 9 . 

c 3 

(1) -Q 

Again we notice that A^ -(y ) = 5.07x10 is very close to 

(1.0) = 5.13X10 -9 , whereas U <:L) | 5 0.4*10~ 2 A. (1) (1.0) 

~ -2 ( 1 ) 

- 0.2x10 (1.0) is practically zero, as it has to be 

analytically . 

It is clear by the continuity of A^^ (u) in the parameter 
y that for u c < y _< 1.05 we have A^ 1 ^ > 0, A^ > 0, and 

A^^ <0, so that the interglacial and the ice-covered earth 
are stable for the entire range of y explored, whereas the 
glacial is unstable for the same range of y. Furthermore the 
ice-covered earth is stable also for smaller y. 

There is one further point of view, which, while illumi- 
nating the significance of the neutral stability of u (x) , 
also argues for our linear stability analysis being sufficient 
for concluding on nonlinear stability or instability of the 
steady-state solutions of (4') corresponding to different 
values of y. This viewpoint has to do with the existence of 
a variational principle for (4 1 ). Indeed 

N (u;y) = 0 

is the Euler-Lagrange equation for the extrema of the functional 

1 

J (u; y) =' | 

0 

here p, r are given by (17b,c) and 


{ 


pu^ +. rG(x 


, u) | 


dx 
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G(x,u) - 

where F(x,u) is defined by (7d). 

Clearly the stable solutions u^(x;l), u^(x;l) corres- 
pond to local minima of J(u;l), whereas u 2 (x;l) is a local 
maximum. As u^(x;p), which is a minimum for ]i > p c , 
coallesces with u 2 (xV y) , which is a maximum for y > ~V- C r sit 

y = y , a saddle point u = u (x) obtains, whereas Uo(x;y ) 

C C J c 

is still a minimum. 

This variational interpretation makes it very plausible 
that solutions u(x,t;y) of the "flow" 


I 


F(x,u) dto 


u t - W (u; y) f U > U c r 

with initial conditions near u^(x;y), that is at a finite 

but small, rather than infinitesimal, distance from Uj(x?y), 

j = 1,2,3, will tend as t -*• 00 towards u^ if u^ is a minimum, 

i.e., j = 1,3, and away from it when u^ is a maximum, i.e., 

j = 2. Similarly, for y - t solutions starting near 

u 3^ x;1J c^ w: *-H s till converge to u^ r but solutions starting 

near u (x) , though they may hover for a long time near u , 
c c 

since t^,t 2 -*■ 00 as y ]x , will eventually move away from it, 
along a negative slope of the saddle, and finally tend towards 
the absolute minimum u^. This seems a rather satisfactory, 
although heuristic, explanation of the results of the time- 
dependent computations of Dwyer and Petersen (1973) and of 
Schneider and Gal-Chen (1973) , which we mentioned alreay in 
Section 5. 
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7 . Concluding Remarks 


We studied the zonal ly- and- vertically averaged energy- 
balance climate model governed by equations (1~3) ; these 
equations are based on simple parameter! zations of albedo, 
greenhouse effect and eddy diffusion of heat in terms of 
yearly averaged sea- level temperature, which is the only 
dependent variable of the model. 

Three positive steady-state solutions of the model, 
symmetric with, respect to the equator, were found by 
accurate numerical computations , and apparently no more 
such solutions exist. These steady states can be identified 
with an interglacial climate, approximating very well the 
one prevailing presently on earth, a glacial climate, and 
a climate during which the earth would be completely ice 
covered. The climates obtained were only slightly changed 
when making small changes in the numerical values of the 
coefficients and when making certain changes in the functional 
form of the model’s equations. However, the bounds on the 
values the albedo can take were essential in order to obtain 
these three climates? also linearizing the outgoing planetary 
radiation resulted in a reduction of the number of solutions. 

We then determined the stability of the time-dependent 
solutions of (1*) under small perturbations about the model's 
steady states. This stability was shown to depend on certain 
properties of a comparison function, which was constructed 
numerically. 
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We found that the interglacial and the "deep freeze" 
climate are stable, and that the glacial climate is unstable. 
This means that the first two can obtain, at least approximate- 
ly, as steady states in a physical system governed by equa- 
tions very similar to (1-3), but that the latter cannot; 
the same is true about these climates as limiting steady 
states for time -dependent numerical solutions of such 
equations . 

We further showed how changes in an important parameter, 
the average intensity of the solar radiation , influence the 
s te ady-s tate s olutions of the mode 1 . The dependence on thi s 
parameter of all steady states was shown to be gradual and 
smooth for increases of up to 5% and decreases of up to 
about 2%. However for a critical value of the parameter, 
equal to 98.15% of its present value, the glacial and inter- 
glacial climates coallesced and they disappeared entirely for 
smaller values of the parameter, leaving the ice-covered earth 
as the only possible stable, steady climate of the model. 

This result is important, as it stresses the difference 
between the stability of a steady state with respect to the 
time evolution of a physical system governed by a given, 
fixed equation, and the stability of a steady state with 
respect to changes in a parameter, which determines the 
behavior of the system. For definiteness, let us call the 
former internal stability and the latter external stability; 
we have shown that the "deep freeze" is stable for our system 
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both internally and externally, that the glacial is unstable 
in both senses, and that the interglacial, or "present 
climate," is internally stable, but externally unstable. 

The limitations of equations (1-3) as a model for the 
description of the long-term behavior of the atmosphere- 
ocean-cryosphere sys*tem, and of energy-balance models 
in general, have been discussed extensively. Because of 
these limitations, we believe that the results above 
should not be taken at face value as statements about the 
climate of our planet. These results, however, seem to 
clarify the physical content and mathematical properties 
of such models. Also, the methods used here could be 
helpful in investigating other models, which will include 
more elaborate and reliable parameterizations of the 
physical phenomena governing climate. We further hope that 
insight gained into the behavior of solutions of a certain 
type of model will advance the formulation of other models, 
and that these will come closer to explaining past changes 
in climate and predicting future changes. 
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Table Captions 

Table 1 . Empirical functions appearing in equation (4) . 

The functions Q, b, z are based on data in Tables 1 
and 2 of Sellers (1973) . The functions C, k^ , k^ 
are based on data provided by Dr. T. Gal-Chen (1974, 
personal communication) , and used in the (SV) model 
of Schneider and Gal-Chen (1973). 

Table 2 . Influence of different modifications in the model's 
equation (4) on the number of steady-state solutions 
and the numerical values of these solutions. The 
existing solutions are identified by the temperature 
at the pole, u^ (0) , j = 1,2,3. In case a solution is 
missing, this is indicated by (-) in the corresponding 
row- and- column location. S stands for the coefficients 
being fitted by cubic splines, B for Bernstein poly- 
nomials. A downward arrow (1) to the right of a comment 
indicates that the equation used in the numerical 
experiments reported in all subsequent rows had the 
feature pointed out in that comment. Otherwise comments 
refer only to the row in which they appear. A left 
arrow, x +■ y, means that the quantity x was replaced 
in the equation by the quantity y. The entries given 
with less than five decimal digits resulted from computa- 
tions with lower precision than indicated in the text. 
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Table 3 . Dependence of the steady states u 1 (x; y) ,u 2 (x; y) on the 
normalized average intensity of the solar radiation, y. 

The columns give u(0), the temperature at the pole, 

£^u, the average temperature, and Au = u(l) - u(0), 
the pole- to-equator temperature difference for u^ and 
for u 2 respectively. The values for y = 0.98151822 
correspond to the bifurcating steady state u = u c (x). 
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Figure Captions 


Figure 1 

Comparison of curve fitting by (i) Bernstein poly- 
nomial approximation, indicated by a dash-dot line, and by 
(ii) cubic spline interpolation, indicated by a solid line. 
Bernstein polynomials are not interpolatory and they are 
variation diminishing, i.e., they have the property of 
smoothing out the data; this results in a rather poor 
approximation. Cubic splines are not variation diminishing 
and they are very good approximants . 

Figure la 

For a very smooth function, like u(x), the two approxi- 
mation procedures yield curves very close to each other. 
Figure lb 

For a function of large total variation, like k(x,u(x)), 
the two procedures yield curves which can differ pointwise by 
as much as 50 percent of the average value. 

Figure 2 

Numerically obtained values of u x (1;Uq) as a function 
of Uq = u(0) . 

Figure 2a 

Comparison of the results for equations (7) , in 
which k = k(x,u), with those for equations (7'), in which 
k = k(x,u(x)); the results for (7) are indicated by a solid 
line, those for (7') by a dash-dot line. 
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Figure 2b 

Results of (7 1 ) for -1330K '<_ u{0) _< 300 K. Notice that 
as tig = u(0) tends towards the ends of the interval. 


u x ( 1 ; u Q ) -> +« 


The 


solution u 4 (x), corresponding to the 


negative root of this curve, u^ = -186 K, does not have a 
physical significance. 

Figure 3 

Comparison of the solutions of (4), indicated by a solid 
line, with those of (4 1 ), indicated by a dash-dot line. The 
circles indicate mesh data for u = u(x) . 

Figure 3a 

Values of the solutions (x) , j = 1,2,3, for (4) and 
for (4 1 ). The respective values for j = 1,3 are practically 
indistinguishable, whereas for j = 2 . a slight difference 
exists between the solution of (4) and that of ( 4 ’ ) . 

Figure 3b 

d. 

Values of the derivatives u.(x), j = 1,2,3, for (4) 

ax j 

and for (4’). The differences are larger than in the function 
values themselves . 

Figure 4 

The first eigenfunctions, w^ 1 ^ (x) , j = 1,2,3 , of 

L j = - h N(u) l u=u .- 

J 3 

Figure 5 

The function , obtained by integrating (7 1 ) 

numerically with \i - \i c and with different values of Uq , 
u Q = u ( 0 > > 120 K. 
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Figure 6 


Dependence of the solutions u_. (x;y) , j = 1,2, on the 
parameter y. The two plots for u^(0), &^Uj are very simi- 
lar? the plot for Au^ = Uj(l)-Uj(0) is rather different, 
although it exhibits the same behavior in the neighborhood 
of the critical point c. The circles indicate the values 
actually computed, for y = y^ , 0.982, 0.985, 0.99, 1.00, 


1.01, 1.02, 1.03, 1.04, 1.05. The letter c distinguishes 
the values of the plotted quantities u^ (0) , ' ^ u j 

corresponding to the bifurcating solution u = u c (x). 


Figure 7 

The bifurcating solution u = u c (x). Notice that the 
ice line, which corresponds to u = 273 K, is at about 45° lat 
for this solution, i.e., an ice cover extending beyond this 
latitude would eventually cover the entire earth. 
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and a completely ice-covered earth* We consider also models 
similar to the main one studied, and determine the number of 
their steady states* All the models have albedo continuously 
varying with latitude and temperature, and entirely diffusive 
horizontal heat transfer* The diffusion is taken to be non- 
linear as well as linear. 

We investigate the stability under small perturbations of 

-4:he -main- -mode l-'s -climates* -A -stability criterion- is derived, 

and its application shows that the "present climate" and the 
"deep freeze" are stable, whereas the model’s glacial is 
unstable. A variational principle is introduced to confirm 
the results of this stability analysis* 

We examine the dependence of the number of steady states 
and of their stability on the average solar radiation* The main 
result is that for a sufficient decrease in solar radiation 
(about 2 percent) the glacial and interglacial solutions 
disappear, leaving the ice-covered earth as the only possible 
climate . 
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